Skip to content
Snippets Groups Projects
Commit 91fd7ece authored by rhaase's avatar rhaase
Browse files

added exercise solution for otsu hypothesis testing

parent 68ac0fa3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
11_Hypothesis_testing/exercise_otsu/Otsu_Method_Comparison.ipynb 0 → 100644
+ 405
− 0
View file @ 91fd7ece
%% Cell type:code id: tags:
``` python
import pandas as pd
# Load data
table_otsu_python = pd.read_csv("otsu_python.csv", delimiter=',')
table_otsu_imagej = pd.read_csv("otsu_imagej.csv", delimiter=',')
# read out columns
measurement_python = table_otsu_python["# Area"];
measurement_imagej = table_otsu_imagej["Area"];
print(measurement_python)
print(measurement_imagej)
```
%% Output
0 199.127197
1 469.055176
2 782.623291
3 990.295410
4 1100.158691
5 1150.970459
6 1173.553467
7 1184.692383
8 1190.338135
9 1185.302734
10 1189.270020
11 1192.169189
12 1206.970215
13 1242.980957
14 1277.008057
15 1309.814453
16 1334.838867
17 1371.765137
18 1375.122070
19 1315.155029
20 1172.790527
21 917.663574
22 587.005615
23 224.304199
24 128.479004
Name: # Area, dtype: float64
0 197.601
1 468.292
2 780.640
3 988.770
4 1099.396
5 1148.682
6 1172.485
7 1183.167
8 1188.507
9 1183.472
10 1187.134
11 1189.117
12 1205.750
13 1238.708
14 1275.330
15 1307.678
16 1331.329
17 1369.019
18 1373.444
19 1313.019
20 1171.265
21 915.527
22 584.717
23 222.931
24 127.411
Name: Area, dtype: float64
%% Cell type:markdown id: tags:
First we determine the mean of both measurementsand compare them
%% Cell type:code id: tags:
``` python
import numpy as np
mean_python = np.mean(measurement_python)
mean_imagej = np.mean(measurement_imagej)
print("Python: " + str(mean_python))
print("ImageJ: " + str(mean_imagej))
print("Difference: " + str(mean_python - mean_imagej))
```
%% Output
Python: 1010.858154296875
ImageJ: 1008.93564
Difference: 1.9225142968749651
%% Cell type:markdown id: tags:
Then we have a look at the histograms for the measurement
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
import numpy as np
def draw_histogram(data):
counts, bins = np.histogram(data, bins=10, range=(0,1400))
plt.hist(bins[:-1], bins, weights=counts)
#plt.axis([0, 10, 0, 4])
plt.show()
draw_histogram(measurement_python)
draw_histogram(measurement_imagej)
```
%% Output
%% Cell type:markdown id: tags:
We should also plot them against each other in a scatter plot to get a first impression on they relate to each other
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
# plot our data
plt.plot(measurement_python, measurement_imagej, "*")
# plot another line which corresponds to identidy
plt.plot([0, 1400], [0, 1400])
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
Now let's have a look at the Bland-Altman plot
%% Cell type:code id: tags:
``` python
# A function for drawing Bland-Altman plots
# source https://stackoverflow.com/questions/16399279/bland-altman-plot-in-python
import matplotlib.pyplot as plt
import numpy as np
def bland_altman_plot(data1, data2, *args, **kwargs):
data1 = np.asarray(data1)
data2 = np.asarray(data2)
mean = np.mean([data1, data2], axis=0)
diff = data1 - data2 # Difference between data1 and data2
md = np.mean(diff) # Mean of the difference
sd = np.std(diff, axis=0) # Standard deviation of the difference
plt.scatter(mean, diff, *args, **kwargs)
plt.axhline(md, color='gray', linestyle='--')
plt.axhline(md + 1.96*sd, color='gray', linestyle='--')
plt.axhline(md - 1.96*sd, color='gray', linestyle='--')
```
%% Cell type:code id: tags:
``` python
# draw a Bland-Altman plot
bland_altman_plot(measurement_python, measurement_imagej)
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
The cofindence interval
%% Cell type:code id: tags:
``` python
data1 = measurement_python
data2 = measurement_imagej
mean = np.mean([data1, data2], axis=0)
diff = data1 - data2 # Difference between data1 and data2
md = np.mean(diff) # Mean of the difference
sd = np.std(diff, axis=0) # Standard deviation of the difference
CI = [md - 2 * sd, md+2*sd]
print(CI)
```
%% Output
[0.3143669053595717, 3.530661688390393]
%% Cell type:code id: tags:
``` python
from statsmodels.stats.weightstats import ttost_ind
pval = ttost_ind(measurement_python, measurement_imagej, low=-10, upp=10)
pval
```
%% Output
(0.47051036223078757,
(0.1097795345763236, 0.4565209578658961, 48.0),
(-0.0743754713943518, 0.47051036223078757, 48.0))
%% Cell type:code id: tags:
``` python
from scipy import stats
def equivalence_ttest_rel(x1, x2, limit):
alpha = 0.05
# print(x1)
# print(x2)
# print(x1 - x2)
# less than test:
t1, p_value1 = stats.ttest_1samp(x1 - x2, limit)
if p_value1 / 2 < alpha and t1 < 0:
print("Reject null-hypothesis: Difference is above limit")
# greater than test:
t2, p_value2 = stats.ttest_1samp(x1 - x2, -limit)
if p_value2 / 2 < alpha and t2 > 0:
print("Reject null-hypothesis: Difference is below -limit")
# return the maximum p_value from both test
# representing the worst case scenario
return [
[t1, p_value1],
[t2, p_value2]
]
tolerance = 10
equivalence_ttest_rel(measurement_python, measurement_imagej, tolerance)
```
%% Output
Reject null-hypothesis: Difference is above limit
Reject null-hypothesis: Difference is below -limit
[[-49.213693923012265, 1.289351642345751e-25],
[72.64029810317597, 1.19906617345889e-29]]
%% Cell type:code id: tags:
``` python
```
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