francescoalemanno

using Statistics: median

function ratio_fit(A,B, pignoleria = 1.0)
    ac = A
    cb = B
    #cb * a = ac
    # 1) divido brutalmente le due quantità
    Q = ac./cb
    # 2) ricavo una stima preliminare del coefficiente 
    # tramite l'uso dello stimatore mediana, che è più robusto della media aritmentica

francescoalemanno

macro 𝔼(v,N,code)
    quote
        M=$(esc(N))
        t = (2-1)/M-1
        z = t/sqrt(1.0-t*t)
        J = exp(-z*z/2.0)/sqrt(2.0π) * (1.0 + z*z)^(3.0/2.0)*2/M
        $(esc(v)) = z
        R=$(esc(code))*J
        for i = 2:M
            t = (2i-1)/M-1

francescoalemanno

γ=4/3
N=6000
y=fill(1.0,N)
z=fill(0.0,N)
x=LinRange(0,10,N)|>collect
dx=x[2]-x[1]
y[2]=1-(dx^2)/6+dx^4/120/(γ-1)
z[2]=-(dx^3)/3 + dx^2/30/(γ-1)
for i =2:(N-1)

francescoalemanno

t=0:0.25:10
n=length(t)
alpha=1.0
mu=zeros(n)
Sigma=zeros(n,n)

for i=1:n
    for j=1:n 

francescoalemanno

using PkgTemplates

t=Template(;julia=v"1.5",plugins=[
    
    GitHubActions(; x86=true),
    Documenter{GitHubActions}(),
])

francescoalemanno

using ApproxBayes
using Distributions

function normaldist(params, constants, targetdata)
  simdata = rand(Normal(params...), 1000)
  ApproxBayes.ksdist(simdata, targetdata), 1
end

using Random
Random.seed!(1)

francescoalemanno

struct Dual{N,T <: Number} <: Number
    x::T
    ϵ::NTuple{N,T}
end

Dual(x,ϵ...)=Dual(x,ϵ)
Dual(x,ϵ::Tuple)=(T=promote(x,ϵ...); Dual(T[1],T[2:end]))
Dual{N,T}(a::Number) where {N,T} = Dual{N,T}(convert(T,a),ntuple(i->zero(T),Val(N)))
Dual{N,T}(a::Dual{N}) where {N,T} = Dual{N,T}(convert(T,a.x),ntuple(i->convert(T,a.ϵ[i]),Val(N)))

francescoalemanno

using ForwardDiff: Dual, value, partials
using PyPlot
using DifferentialEquations

f(β,m)=tanh(β*m)-m
g(β,m,β0,m0)=f(β,m)-f(β0,m0)

function mdot(β,m,β0,m0)
    dm=Dual{1}(m,oftype(m,1))
    dβ=Dual{1}(β,oftype(β,1))

francescoalemanno

using ForwardDiff: Dual, value, partials, npartials

function test(μ,σ,iters)
    x=(μ+tanh(randn())*σ)^2
    N=1
    xsq=x*x
    for i in 1:iters
        v=(μ+tanh(randn())*σ)^2
        x+=v
        xsq+=v*v

francescoalemanno

using PyCall
using PyPlot

function matplotlibstyle()
    py"""
    from matplotlib import rcParams
    import numpy as np
    def mplfigsize(hscale, publicationcols=1,columnsinches=3.1):
        bonuscol=-0.12
        W=columnsinches*(publicationcols+bonuscol)