

In [1]:

    
import PyPlot.plt
using NLopt









    



INFO: Loading help data...

Optimization

This notebook is designed to figure out if we can use NLopt as a first pass to fit the disks.

To figure out how to use it, let's try doing a simple Gaussian example



In [2]:

    
n = 100









    Out[2]:





100



In [13]:

    
xx = linspace(0, 10., n) #vector of x points;



In [4]:

    
# Compute a Gaussian with amplitude, mean, and variance
# p = [amp, mean, variance] 
function Gauss(x::Vector{Float64}, p::Vector{Float64})
    p[1] * exp(-0.5 * (x - p[2]).^2 ./ p[3]^2)    
end









    Out[4]:





Gauss (generic function with 1 method)



In [5]:

    
truth = Float64[2., 4.0, 0.8]









    Out[5]:





3-element Array{Float64,1}:
 2.0
 4.0
 0.8



In [14]:

    
yy = Gauss(xx, truth) + randn(n) #add some noise to the Gaussian;



In [15]:

    
plt.plot(xx, yy, "o")
plt.show()



In [19]:

    
#Define the objective function that we want to minimize
function chi2(p::Vector, grad::Vector)
    println("parameters ", p)
    sumabs2(yy - Gauss(xx, p))
end

function chi2(p::Vector)
    println("parameters ", p)
    sumabs2(yy - Gauss(xx, p))
end









    Out[19]:





chi2 (generic function with 2 methods)



In [16]:

    
chi2(truth, zeros(2))









    



parameters [2.0,4.0,0.8]






    Out[16]:





104.55200238319387



In [10]:

    
opt = Opt(:LN_COBYLA, 3)









    Out[10]:





Opt(:LN_COBYLA, 3)



In [12]:

    
min_objective!(opt, chi2)
xtol_rel!(opt,1e-4)



In [17]:

    
(optf,optx,ret) = optimize(opt, [1.0, 2.0, 1.0])









    



parameters [1.0,2.0,1.0]
parameters [2.0,2.0,1.0]
parameters [1.0,4.0,1.0]
parameters [1.0,4.0,2.0]
parameters [0.42944465691686184,5.599377267204626,2.1869884749564754]
parameters [1.4823434147977959,4.256481308708498,2.0299860763021007]
parameters [0.9497084687237245,3.0132213835229056,2.0635433111961863]
parameters [0.7779413841985976,4.1832801252309775,2.069224779322107]
parameters [0.9819244050365715,3.963487454974981,1.7513235407589471]
parameters [1.103381162094574,3.976788675725573,1.722530737029727]
parameters [1.174297914250884,3.773310291872704,1.7068798236456244]
parameters [1.2798228781073848,3.655495933880015,1.6749513682213983]
parameters [1.2556790191319964,3.439562307557534,1.6167685414631325]
parameters [1.2719861318641534,3.649043884003809,1.6841481949246835]
parameters [1.4007634982884403,3.59586950410819,1.6854174376071844]
parameters [1.2870496584680249,3.699819040344671,1.6612473574122553]
parameters [1.2188676779159382,3.8410267844622,1.5838437994946475]
parameters [1.28825752624537,3.7248217219619044,1.6892117894249472]
parameters [1.2235909007259622,3.4848434259470857,1.6545872245735338]
parameters [1.2556492361652525,3.7738899004628017,1.776431516043326]
parameters [1.2927074860995822,3.7315140076662927,1.6010551963934654]
parameters [1.2906015379766078,3.803578491257457,1.5500307656402659]
parameters [1.3201564185434556,3.745094366423773,1.6049694125611804]
parameters [1.3381255189440178,3.6451664111269637,1.5720003002233538]
parameters [1.3571129376912139,3.542643305846831,1.5417020964771193]
parameters [1.4194878319102016,3.5467967655578185,1.5450650957749763]
parameters [1.3568862015316374,3.5182130047456854,1.5639873733741798]
parameters [1.3570651060132253,3.5617812098242005,1.479938977466104]
parameters [1.3539465128757335,3.6759385491014958,1.4546694335619312]
parameters [1.3535542282865887,3.636165296032747,1.3954189635394627]
parameters [1.2978102301318977,3.6235025478413063,1.3678731767147985]
parameters [1.4145242692861608,3.6383651417560943,1.381718807693311]
parameters [1.4708729424672788,3.6365806767845403,1.3546950967727345]
parameters [1.4727381732324005,3.5983591258643273,1.2952177971071794]
parameters [1.4687301084706426,3.6756996041183365,1.3482442664643173]
parameters [1.5213193399683458,3.5727043641255873,1.3731711664153363]
parameters [1.482195524495028,3.6933395607836177,1.3481400294184116]
parameters [1.4685701423895121,3.622379440624758,1.340968518186453]
parameters [1.4671550451336923,3.623175865663474,1.3409389519151502]
parameters [1.4727866278889847,3.5925955554501128,1.3431113949931364]
parameters [1.476133107987496,3.5936299200701103,1.3409024732267458]
parameters [1.486483330384611,3.5936072528179412,1.3291972012145634]
parameters [1.5009116961102071,3.6024211940566713,1.333264325192081]
parameters [1.514921537234759,3.600226900690047,1.3264334214583051]
parameters [1.519898412643348,3.5784163750255753,1.3364557255211063]
parameters [1.5155710728796867,3.630844356389394,1.323373751529931]
parameters [1.5117439990369057,3.597500934206707,1.3194276606906148]
parameters [1.510953893758055,3.5928977739572443,1.3120038059757262]
parameters [1.5111049521067446,3.597104215540497,1.3198423478569727]
parameters [1.50956665665716,3.6122775337833897,1.3181209998130505]
parameters [1.5127891218329015,3.6263323558041063,1.319246330401201]
parameters [1.5102032382772386,3.6122202698792765,1.316303064911454]
parameters [1.5032189953263333,3.6143199507986274,1.3129634030031228]
parameters [1.5144821418982812,3.6119558564277807,1.3097678784999203]
parameters [1.5158763970248992,3.604910680536279,1.310719832211427]
parameters [1.5157898840844535,3.612010000738389,1.3104930197033626]
parameters [1.5169431689432742,3.6193698712974456,1.311114951721672]
parameters [1.5144575373127245,3.6129094904747143,1.3091375570889576]
parameters [1.516080521570233,3.61198838250364,1.3099293730158459]
parameters [1.516452395194084,3.6152686997688632,1.3113421348030034]
parameters [1.5167548892515914,3.6146635265099443,1.3122199804575756]
parameters [1.516032308256197,3.615401152900353,1.3127636051889338]
parameters [1.5160520345547808,3.615930665817396,1.312798094160589]
parameters [1.5155589023613905,3.6136929100588837,1.3127697512291978]
parameters [1.5161687052463881,3.6122800407771662,1.31305744004728]
parameters [1.5156503162849408,3.6136469120650223,1.3125872216170296]
parameters [1.514675716217803,3.6132234020887837,1.3124254445073862]
parameters [1.5154669400334069,3.61371170811045,1.313189835312352]
parameters [1.5154013884625677,3.6155772912856237,1.3125669199294399]
parameters [1.515714696473193,3.6127755696275083,1.3127083638765078]
parameters [1.5155367806037254,3.614150722531741,1.3126877926571185]
parameters [1.5156454378953474,3.613712762096418,1.312782673586152]
parameters [1.5153728822533405,3.613465872608586,1.3126596845983913]
parameters [1.5151740889442644,3.6134347823692274,1.3128005540922363]
parameters [1.5150278497713456,3.613776776912285,1.3127057981836647]
parameters [1.514801842468558,3.613613093015097,1.3126630519638038]
parameters [1.5149586140996578,3.614091168720106,1.3128792931816764]
parameters [1.515043094501242,3.6138598844090786,1.3126161272061718]
parameters [1.5150302930659294,3.613738563188123,1.3126967902631692]
parameters [1.5150436877832518,3.613650467063999,1.3127855615973736]
parameters [1.514948274555952,3.613604567028346,1.3128047869124029]
parameters [1.514878791832542,3.6135183038253094,1.3128623352929012]
parameters [1.5148781230532944,3.613354318494234,1.3128050923410874]
parameters [1.5148267141310823,3.613684384115677,1.3128821405270206]
parameters [1.5148612403392567,3.6137525268327555,1.3129695881311982]






    Out[17]:





(96.58404355224967,[1.51486,3.61375,1.31297],:XTOL_REACHED)

The MCMC approach



In [18]:

    
using MCMC



In [20]:

    
mcmodel = model(chi2, init=[1.0, 2.0, 1.0])









    



parameters [1.0,2.0,1.0]






    Out[20]:





LikelihoodModel, with 1 parameter(s)



In [21]:

    
mcmodel.eval(mcmodel.init)









    



parameters [1.0,2.0,1.0]






    Out[21]:





129.06597228854199



In [27]:

    
mcrunner = SerialMC(nsteps=100, burnin=10)









    Out[27]:





SerialMCBaseRunner(10,1,100,11:1:100,false)



In [41]:

    
mcsampler = MH()









    Out[41]:





MH(true,nothing,(anonymous function))



In [32]:

    
sigma = [0.2, 0.2, 0.1]









    Out[32]:





3-element Array{Float64,1}:
 0.2
 0.2
 0.1



In [42]:

    
mcchain = run(mcmodel, mcsampler, mcrunner)









    



parameters [1.0,2.0,1.0]






    



`convert` has no method matching convert(::Type{ScalMat}, ::Float64)
while loading In[42], in expression starting on line 1

 in schedule_and_wait at task.jl:256
 in consume at task.jl:167
 in anonymous at no file:26
 in run at /home/ian/.julia/MCMC/src/runners/SerialMC.jl:40
 in run at /home/ian/.julia/MCMC/src/runners/SerialMC.jl:29



In [ ]: