Inverse function in modeling

Modeling
10

Great! I implemented models by build-in Powell’s hybrid algorithm and Newton’s method to check the speed, . The latter model is extremely fast!! This solves my problem completely. I attach the scripts for someone who are interested in.

I can’t thank you enough!

Stan

Build-in function

model = """
functions{
    // system of inverse Batschelet transformation
    vector invBatschelet(vector y_init, vector t, real[] x_r, int[] x_i){
        // Get length of array
        int I = x_i[1];
        // value of interest
        vector[I] y ;
        // parameters (take it as fixed value would be understandable)
        vector[I] theta = t[:I] ;
        real nu = t[I+1] ;
        
        // construct equation to be zero
        y = -theta + y_init - nu * (1 + cos(y_init)) ;
        return y ;
    }
}

data {
    int I ;
    vector<lower=-pi(), upper=pi()>[I] RADIAN ;
    int K ;
}

transformed data {
    real x_r[0] ;
    int x_i[1] ;
    vector<lower=-pi(), upper=pi()>[I] THETA_INIT ;
    vector<lower=0.0>[K] A ;
    
    // set integer parameter to pass solver
    x_i[1] = I ;
    // set initial value of inverse transformed value
    for (i in 1:I){
        THETA_INIT[i] = 0 ;
    }
    for (k in 1:K){
        A[k] = 50 / K ;
    }
}

parameters {
    simplex[K] alpha ;
    unit_vector[2] O ;
    real<lower=0.0> kappa[K] ;
    real<lower=-1.0, upper=1.0> nu[K] ;
}

transformed parameters{
    real mu ;
    
    mu = atan2(O[1], O[2]) ;
}

model {
    alpha ~ dirichlet(A) ;
    for(k in 1:K){
        vector[I] THETA ;
        vector[I+1] t ;
        for (i in 1:I){
            // transform angle by location parameter
            t[i] = RADIAN[i] - mu ;
            if (t[i] > pi()){
                t[i] = t[i] - 2*pi() ;
            }else if (t[i] <= -pi()){
                t[i] = t[i] + 2*pi() ;
            }
        }
        t[I+1] = nu[k] ;
        THETA = algebra_solver(invBatschelet, THETA_INIT, t, x_r, x_i) ;   
        THETA ~ von_mises(0, kappa[k]) ;
    }
}
"""

model = StanModel(model_code=model)

Newton’s method

model = """
functions{
    vector inv_transformation_Batschelet(vector theta, real mu, real nu, int I){
        vector[I] t;
        t = theta ;
        // Number of iteration was evaluated with scipy's simulation
        for (i in 1:8){
            t = t - ((t - nu*(1+cos(t-mu)) - theta) ./ (1 + nu * (sin(t-mu)))) ;
        }
        return t;
    }
}

data {
    int I ;
    vector<lower=-pi(), upper=pi()>[I] RADIAN ;
    int K ;
}

transformed data {
    vector<lower=0.0>[K] A ;

    for (k in 1:K){
        A[k] = 50 / K ;
    }
}

parameters {
    simplex[K] alpha ;
    unit_vector[2] O ;
    real<lower=0.0> kappa[K] ;
    real<lower=-1.0, upper=1.0> nu[K] ;
}

transformed parameters{
    real mu ;
    
    mu = atan2(O[1], O[2]) ;
}

model {
    alpha ~ dirichlet(A) ;
    for(k in 1:K){
        vector[I] THETA ;
        THETA = inv_transformation_Batschelet(RADIAN, mu, nu[k], I) ;
        THETA ~ von_mises(mu, kappa[k]) ;
    }
}
"""

model = StanModel(model_code=model)

Python

Code

# Change sample size with your demand
I = 3000
RADIAN = vonmises.rvs(size=I, loc=0, kappa=0.3)
stan_data = {
    "I": I,
    "RADIAN": RADIAN,
    "K": 1
}

%timeit fit = model.optimizing(stan_data, seed=1234, init_alpha=1e-10)

Result

time_methods