Last active
November 29, 2018 23:01
-
-
Save ellielinc/fce4bb2e14d3653ede5422b50666f441 to your computer and use it in GitHub Desktop.
Uses multiple gridsearches to minimize the chi square of the best-fit diameter and limb darkening coefficient of a star with the linear visibililty function and generates errorbars for the best-fit parameters by calculating the probability density of those parameters.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
function errorbar_LDlin_prob(data::OIdata, | |
dataname="", diam_min=1.0, diam_max=9.0, LD_min=0.01, LD_max=0.9, iter=200) | |
chisq=param->norm((abs2.(visibility_ldlin(param,data.v2_baseline))-data.v2)./data.v2_err).^2; | |
###taken from bic.jl in numericalmethods | |
function logsumexp(X) | |
# Compute log(sum(exp(X))) while avoiding numerical over/underflow. | |
return log(sum(exp.(X-maximum(X))))+maximum(X) | |
end | |
#initial grid search | |
LD = linspace(LD_min,LD_max,iter); | |
diam = linspace(diam_min, diam_max, iter); | |
chi2val1 = zeros(length(diam),length(LD)); | |
for i=1:length(diam) | |
for j=1:length(LD) | |
chi2val1[i,j]=chisq([diam[i],LD[j]]); | |
end | |
end | |
#chi2 map | |
indx = find(chi2val1.>(minimum(chi2val1)+1000)); | |
chi2val_filtered = copy(chi2val1); | |
chi2val_filtered[indx] = 0; | |
minchisq1 = find(chi2val_filtered.>0); | |
chi2_1 = minimum(chi2val_filtered[minchisq1]); | |
fig1 = imshow(rotl90(chi2val_filtered./data.nv2),ColorMap("magma"), extent=[diam[1], diam[end], 0, 1]) | |
colorbar(); title("Initial Chi2 Map $dataname"); | |
println("Best Chi2 Value of initial $dataname gridsearch = $chi2_1") | |
#refine params | |
best1=ind2sub(chi2val1,indmin(chi2val1)); | |
diam_min=diam[best1[1]]-1; | |
diam_plus=diam[best1[1]]+1; | |
LD_min=LD[best1[2]]-0.1 | |
LD_plus=LD[best1[2]]+0.1 | |
#2nd gridsearch | |
LD = linspace(LD_min, LD_plus, iter); | |
diam = linspace(diam_min, diam_plus, iter); | |
chi2val2=zeros(length(diam),length(LD)); | |
for i=1:length(diam) | |
for j=1:length(LD) | |
chi2val2[i,j]=chisq([diam[i],LD[j]]); | |
end | |
end | |
#refined chi2 map | |
indx = find(chi2val2.>(minimum(chi2val2)+500)); | |
chi2val_filtered = copy(chi2val2); | |
chi2val_filtered[indx] = 0; | |
minchisq2 = find(chi2val_filtered.>0) | |
chi2_2 = minimum(chi2val_filtered[minchisq2]) | |
figure() | |
fig2 = imshow(rotl90(chi2val_filtered./data.nv2),ColorMap("magma"), extent=[diam[1], diam[end], 0, 1]) | |
colorbar(); title("Refined Chi2 Map $dataname"); | |
println("Best Chi2 Value of refined $dataname gridsearch = $chi2_2") | |
#best diam param | |
best = ind2sub(chi2val2,indmin(chi2val2)); | |
d_slice = diam[find(mapslices(logsumexp, -chi2val_filtered./2, dims=[2]).>0)] | |
diam_best = diam[best[1]] | |
diam_best_min = diam_best-minimum(d_slice) | |
diam_best_max = maximum(d_slice)-diam_best | |
LD_slice = LD[find(mapslices(logsumexp, -chi2val_filtered./2, dims=[1]).>0)] | |
LD_best = LD[best[2]] | |
LD_best_min = LD_best-minimum(LD_slice) | |
LD_best_max = maximum(LD_slice)-LD_best | |
bestchi2 = minimum(chi2val2)./data.nv2 | |
println("Best Chi2: $bestchi2") | |
println("Best Diameter: $diam_best, +$diam_best_max, -$diam_best_min") | |
println("Best LDC: $LD_best, +$LD_best_max, -$LD_best_min") | |
return bestchi2, diam_best, diam_best_max, diam_best_min, LD_best, LD_best_max, LD_best_min | |
end |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment