Skip to content
Snippets Groups Projects
Commit 873e5a29 authored by Philipp Frank's avatar Philipp Frank
Browse files

Merge branch 'iid2022' into 'main' 

make evidences an array

See merge request !5
parents 86057d21 c720aeae
No related branches found
No related tags found
1 merge request!5make evidences an array
Pipeline #148616 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
demo_poisson.ipynb
+ 4
− 4
View file @ 873e5a29
%% Cell type:markdown id: tags:
Nifty tutorial for Poisson count data
=====================================
%% Cell type:markdown id: tags:
Setup
-----
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import nifty8 as ift
from utils import plot_2D, load_psf, geovi_sampling, plot_posterior
ift.random.push_sseq_from_seed(42)
evidences = []
evidences = np.zeros(3)
```
%% Cell type:code id: tags:
``` python
# Load data and visualize
data = np.load('data/poisson.npz')
```
%% Cell type:code id: tags:
``` python
position_space = ift.RGSpace([128, 128])
# Homogeneous poisson process
print(model1)
```
%% Cell type:code id: tags:
``` python
# Set up likelihood & PSF
```
%% Cell type:code id: tags:
``` python
# Inference model 1
print(evidence)
evidences += [evidence, ]
evidences[0] = evidence
```
%% Cell type:code id: tags:
``` python
# Independent poisson process
# Inference model 2
evidences += [evidence, ]
evidences[1] = evidence
```
%% Cell type:code id: tags:
``` python
# Compare evidence
print(evidences)
```
%% Cell type:code id: tags:
``` python
# Diffuse poisson process
args = {
'offset_mean': .5,
'offset_std': (1., 1E-5),
# Amplitude of field fluctuations
'fluctuations': (1.5, 0.5), # 1.0, 1e-2
# Exponent of power law power spectrum component
'loglogavgslope': (-4., 1), # -6.0, 1
# Amplitude of integrated Wiener process power spectrum component
'flexibility': (1., 0.2), # 2.0, 1.0
# How ragged the integrated Wiener process component is
'asperity': (0.1, 0.01), # 0.1, 0.5
# Name of the input keys
'prefix' : 'diffuse'
}
correlated_field = ift.SimpleCorrelatedField(position_space, **args)
pspec = correlated_field.power_spectrum
```
%% Cell type:code id: tags:
``` python
# Prior samples
```
%% Cell type:code id: tags:
``` python
# Inference model 3
evidences += [evidence, ]
evidences[2] = evidence
```
%% Cell type:code id: tags:
``` python
# Compare evidence
print(evidences)
print(evidences - evidences[-1])
```
%% Cell type:markdown id: tags:
Posterior visualization
-----------------------
%% Cell type:code id: tags:
``` python
plot_posterior(samples, data, model3, diffuse, model2, pspec)
```
......
demo_poisson_solution.ipynb
+ 5
− 5
View file @ 873e5a29
%% Cell type:markdown id: tags:
Nifty tutorial for Poisson count data
=====================================
%% Cell type:markdown id: tags:
Setup
-----
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import nifty8 as ift
from utils import plot_2D, load_psf, geovi_sampling, plot_posterior
ift.random.push_sseq_from_seed(42)
evidences = []
evidences = np.zeros(3)
```
%% Cell type:code id: tags:
``` python
# Load data and visualize
data = np.load('data/poisson.npz')
print(data['data'].shape)
plot_2D(data['data'], 'Data')
```
%% Cell type:code id: tags:
``` python
position_space = ift.RGSpace([128, 128])
# Homogeneous poisson process
projection = ift.VdotOperator(ift.full(position_space, 1.)).adjoint
model1 = ift.FieldAdapter(projection.domain, 'hom')
model1 = ift.exp(5. * model1)
model1 = projection @ model1
print(model1)
```
%% Cell type:code id: tags:
``` python
# Set up likelihood & PSF
d = ift.makeField(position_space, data['data'])
likelihood = ift.PoissonianEnergy(d)
PSF_op, psf = load_psf(position_space)
plot_2D(psf, 'PSF')
likelihood = likelihood @ PSF_op
```
%% Cell type:code id: tags:
``` python
# Inference model 1
samples, evidence = geovi_sampling(likelihood @ model1)
plot_2D(samples.average(model1).val, 'model1 posterior mean')
print(evidence)
evidences += [evidence, ]
evidences[0] = evidence
```
%% Cell type:code id: tags:
``` python
# Independent poisson process
model2 = ift.InverseGammaOperator(position_space, 2., 3.).ducktape('independent')
# Inference model 2
samples, evidence = geovi_sampling(likelihood @ model2)
plot_2D(samples.average(model2).val, 'model2 posterior mean')
print(evidence)
evidences += [evidence, ]
evidences[1] = evidence
```
%% Cell type:code id: tags:
``` python
# Compare evidence
print(evidences)
print(evidences[:2])
```
%% Cell type:code id: tags:
``` python
# Diffuse poisson process
args = {
'offset_mean': .5,
'offset_std': (1., 1E-5),
# Amplitude of field fluctuations
'fluctuations': (1.5, 0.5), # 1.0, 1e-2
# Exponent of power law power spectrum component
'loglogavgslope': (-4., 1), # -6.0, 1
# Amplitude of integrated Wiener process power spectrum component
'flexibility': (1., 0.2), # 2.0, 1.0
# How ragged the integrated Wiener process component is
'asperity': (0.1, 0.01), # 0.1, 0.5
# Name of the input keys
'prefix' : 'diffuse'
}
correlated_field = ift.SimpleCorrelatedField(position_space, **args)
pspec = correlated_field.power_spectrum
diffuse = correlated_field.exp()
model3 = diffuse + model2
```
%% Cell type:code id: tags:
``` python
# Prior samples
pl = ift.Plot()
for _ in range(9):
pl.add(model3(ift.from_random(model3.domain)))
pl.output()
```
%% Cell type:code id: tags:
``` python
# Inference model 3
samples, evidence = geovi_sampling(likelihood @ model3)
plot_2D(samples.average(model3).val, 'model3 posterior mean')
print(evidence)
evidences += [evidence, ]
evidences[2] = evidence
```
%% Cell type:code id: tags:
``` python
# Compare evidence
print(evidences)
print(evidences - evidences[-1])
```
%% Cell type:markdown id: tags:
Posterior visualization
-----------------------
%% Cell type:code id: tags:
``` python
plot_posterior(samples, data, model3, diffuse, model2, pspec)
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment