1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
#!/usr/bin/env python3
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""
Example of a script that performs histogram analysis of an activation image.
This is based on a real fMRI image.
Simply modify the input image path to make it work on your preferred image.
Needs matplotlib
Author : Bertrand Thirion, 2008-2009
"""
import os
import numpy as np
import scipy.stats as st
try:
import matplotlib.pyplot as plt
except ImportError:
raise RuntimeError("This script needs the matplotlib library")
# Local import
from get_data_light import DATA_DIR, get_second_level_dataset
from nibabel import load
import nipy.algorithms.statistics.empirical_pvalue as en
# parameters
verbose = 1
theta = float(st.t.isf(0.01, 100))
# paths
mask_image = os.path.join(DATA_DIR, 'mask.nii.gz')
input_image = os.path.join(DATA_DIR, 'spmT_0029.nii.gz')
if (not os.path.exists(mask_image)) or (not os.path.exists(input_image)):
get_second_level_dataset()
# Read the mask
nim = load(mask_image)
mask = nim.get_fdata()
# read the functional image
rbeta = load(input_image)
beta = rbeta.get_fdata()
beta = beta[mask > 0]
mf = plt.figure(figsize=(13, 5))
a1 = plt.subplot(1, 3, 1)
a2 = plt.subplot(1, 3, 2)
a3 = plt.subplot(1, 3, 3)
# fit beta's histogram with a Gamma-Gaussian mixture
bfm = np.array([2.5, 3.0, 3.5, 4.0, 4.5])
bfp = en.gamma_gaussian_fit(beta, bfm, verbose=1, mpaxes=a1)
# fit beta's histogram with a mixture of Gaussians
alpha = 0.01
pstrength = 100
bfq = en.three_classes_GMM_fit(beta, bfm, alpha, pstrength,
verbose=1, mpaxes=a2)
# fit the null mode of beta with the robust method
efdr = en.NormalEmpiricalNull(beta)
efdr.learn()
efdr.plot(bar=0, mpaxes=a3)
a1.set_title('Fit of the density with \n a Gamma-Gaussian mixture')
a2.set_title('Fit of the density with \n a mixture of Gaussians')
a3.set_title('Robust fit of the density \n with a single Gaussian')
plt.show()
|
