Fitting stream with varying RAR (theta_0, W_0) and disk parameters.

Dear @Charly-Arguelles @danielcarpintero @ankrut :

I quote from the previous issue:

In my opinion, one remaining point would bee to check (at least minimally) what can be gained further regarding the fitting power of the rotation curve, by varying somewhat the baryons (mainly the disk parameters). But of course, it should be re-done the analysis for GD-1, in a kind of iterative fashion (because we are not doing all the fittings in parallel as we have already agreed in the past). Martin, one question from my side (which I forgot): are we using the Miyamoto Nagai disk right?, and, what about the baryonic disks (formula and specific free parameters) used by the other authors we are comparing with in the Vrot curve?. To me, there is certainly some extra freedom about the baryons yet, in the sense that, for example the contribution of the disk could make the RAR DM tail start to decrease further in radius..

I will do some fits varying the barions as follows: 15% for central values of the masses and 10% for characteristic scales, according to the first paragraph of section 3.1.2 of this paper

@ankrut @Charly-Arguelles @danielcarpintero , here the results.

The varying parameters in the fit are theta_0, d_theta (=W_0-theta_0) plus the thin and thick disk parameters (Please note that I use the same factors for both disks). The disk parameters are the total mass, the r-scale and the z-scale.

def pot_model(ener_f, theta_0, W_0, beta_0, alpha, scale, scaleb):
    """Potential model including barions and dark matter halo."""
    # Set barionic potentials
    M_gal = 2.32e7  # Msun
    bulge = pot.Plummer(460.0*M_gal, 0.3)
    thin = pot.MiyamotoNagai(1700.0*M_gal*alpha, 5.3*scale, 0.25*scaleb)
    thick = pot.MiyamotoNagai(1700.0*M_gal*alpha, 2.6*scale, 0.8*scaleb)
    # Set halo potential
    param = np.array([ener_f, theta_0, W_0, beta_0])
    halo = pot.RAR(param)
    return bulge, thin, thick, halo

The a priori bounds for those parameters (theta_0, d_theta, alpha, scale, scaleb) are: # bounds = ((34, 39), (25, 30), (0.85, 1.15), (0.9, 1.1), (0.9, 1.1)) I made two runs with a rather large number of iterations-population equal to 40-40 and 100-100, respectively. In the first case I obtained: theta_0, d_theta, alpha, scale, scaleb = 3.642214195068059013e+01, 2.758970214272743249e+01 1.089703158216322576e+00 1.033675113036105708e+00 9.062517097351714401e-01

In the second case I obtained: 3.626724510339099083e+01 2.747475602895485736e+01 1.148997022079204511e+00 1.063196394644405185e+00 9.072897756291734561e-01

Both cases are compatible between them and also in agreement with the case with fixed barions (alpha=1, scale=1, scaleb=1). The chi^2 reduces a bit to 27 in this two cases with varying disk. In the case of fixed disk, the chi^2 is 29. This slight difference is not noticeable from the stream plots. Please notice that in this case I am using beta_0=1.25e-5 while in the original fits I used 1.2e-5. Here is a plot of the rotation curve for the second case, rotation_curve_beta0_1 25e-5 in order to compare with the plots in the previous issue.

How do you compute chi^2 in these experiments?

El mié, 2 mar 2022 a las 16:04, martinmestre @.***>) escribió:

@ankrut https://github.com/ankrut @Charly-Arguelles https://github.com/Charly-Arguelles @danielcarpintero https://github.com/danielcarpintero , here the results.

The varying parameters in the fit are theta_0, d_theta (=W_0-theta_0) plus the thin and thick disk parameters (Please note that I use the same factors for both disks). The disk parameters are the total mass, the r-scale and the z-scale.

def pot_model(ener_f, theta_0, W_0, beta_0, alpha, scale, scaleb): """Potential model including barions and dark matter halo."""

Set barionic potentials
M_gal = 2.32e7  # Msun
bulge = pot.Plummer(460.0*M_gal, 0.3)
thin = pot.MiyamotoNagai(1700.0*M_gal*alpha, 5.3*scale, 0.25*scaleb)
thick = pot.MiyamotoNagai(1700.0*M_gal*alpha, 2.6*scale, 0.8*scaleb)
# Set halo potential
param = np.array([ener_f, theta_0, W_0, beta_0])
halo = pot.RAR(param)
return bulge, thin, thick, halo
The a priori bounds for those parameters (theta_0, d_theta, alpha, scale, scaleb) are:

bounds = ((34, 39), (25, 30), (0.85, 1.15), (0.9, 1.1), (0.9, 1.1))

I made two runs with a rather large number of iterations-population equal to 40-40 and 100-100, respectively. In the first case I obtained: theta_0, d_theta, alpha, scale, scaleb = 3.642214195068059013e+01, 2.758970214272743249e+01 1.089703158216322576e+00 1.033675113036105708e+00 9.062517097351714401e-01

In the second case I obtained: 3.626724510339099083e+01 2.747475602895485736e+01 1.148997022079204511e+00 1.063196394644405185e+00 9.072897756291734561e-01

Both cases are compatible between them and also in agreement with the case with fixed barions (alpha=1, scale=1, scaleb=1). The chi^2 reduces a bit to 27 in this two cases with varying disk. In the case of fixed disk, the chi^2 is 29. This slight difference is not noticeable from the plots. Please notice that in this case I am using beta_0=1.25e-5 while in the original fits I used 1.2e-5. Here is a plot of the rotation curve for the second case, [image: rotation_curve_beta0_1 25e-5] https://user-images.githubusercontent.com/58749264/156430375-774527c5-8a6f-496a-a89f-8a00786af3d2.png in order to compare with the plots in the previous issue.

— Reply to this email directly, view it on GitHub https://github.com/martinmestre/stream-fit/issues/5#issuecomment-1057277705, or unsubscribe https://github.com/notifications/unsubscribe-auth/ATTASC6DMI3WHO64QWBWIBTU563SJANCNFSM5PYJNEUQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you were mentioned.Message ID: @.***>

Excellent @martinmestre

@ankrut , @danielcarpintero , as we were discussing today with Martin (from what he figured it out), by defining a NEW likelihood, where we add the \Chi² with the Vrot data (e.g. Sofue+ 2020), to the already used \Chi² of the stream data, it would be very interesting to see the new result of maximization of this new Likelihood, because now we may obtain a considerably better fit to the V_rot , while still providing a good fit to the stream.

As Daniel asks, it would be nice if you can write down how you are planning to define the Lieklihhods (either in the original case, as in the new one) Carlos

@danielcarpintero :

How do you compute chi^2 in these experiments?

Like this:

def chi2(w_0, ener_f, beta_0, ic):
    """Chi^2 function."""
    import wrap
    theta_0 = w_0[0]
    d_theta = w_0[1]
    alpha = w_0[2]
    scale = w_0[3]
    scaleb = w_0[4]
    W_0 = theta_0 + d_theta
    pot_list = pot_model(ener_f, theta_0, W_0, beta_0, alpha, scale, scaleb)
    phi_1, phi_2, d_hel, mu_ra, mu_dec, v_hel, x, y, z, v_circ = orbit_model(
        ic[0], ic[1], ic[2], ic[3], ic[4], ic[5], pot_list)
    cfg.phi_2_spl = interp1d(phi_1, phi_2, kind='cubic')
    cfg.d_hel_spl = interp1d(phi_1, d_hel, kind='cubic')
    cfg.v_hel_spl = interp1d(phi_1, v_hel,  kind='cubic')
    cfg.mu_ra_spl = interp1d(phi_1, mu_ra, kind='cubic')
    cfg.mu_dec_spl = interp1d(phi_1, mu_dec, kind='cubic')
    cfg.phi_1_min = np.amin(phi_1.value)
    cfg.phi_1_max = np.amax(phi_1.value)

    sum = np.zeros(5)

    y_mod = wrap.phi_2_wrap(Iba_sky['phi_1'])
    y_dat = Iba_sky['phi_2']
    sigma2 = 0.5**2
    sum[0] = np.sum((y_dat-y_mod)**2 / sigma2)

    y_mod = wrap.d_hel_wrap(Iba_sky['phi_1'])
    y_dat = Iba_sky['d_hel']
    sigma2 = 1.5**2
    sum[1] = np.sum((y_dat-y_mod)**2 / sigma2)

    y_mod = wrap.mu_ra_wrap(Iba_sky['phi_1'])
    y_dat = Iba_sky['mu_ra']
    sigma2 = 4.0
    sum[2] = np.sum((y_dat-y_mod)**2 / sigma2)

    y_mod = wrap.mu_dec_wrap(Iba_sky['phi_1'])
    y_dat = Iba_sky['mu_dec']
    sigma2 = 4.0
    sum[3] = np.sum((y_dat-y_mod)**2 / sigma2)

    y_mod = wrap.v_hel_wrap(Iba_sky['phi_1'])
    y_dat = Iba_sky['v_hel']
    sigma2 = 100.0
    sum[4] = np.sum((y_dat-y_mod)**2 / sigma2)

    print('chi^2 =', np.sum(sum), ' v_circ = ', v_circ, ' theta_0=', theta_0, ' W_0=', W_0,
          ' alpha=', alpha, ' scale=', scale, 'scaleb=', scaleb)
    return np.sum(sum)

It is a discrete sum of the differences between the Ibata polynomials and our orbits, for the five observables: position in sky (only Phi_2 because Phi_1 is the independent variable here), distance, proper motion and radia velocity, all added together. Likke an L_2 norm between functions.

martinmestre / stream-fit

Fitting stream with varying RAR (theta_0, W_0) and disk parameters. #5

Set barionic potentials

bounds = ((34, 39), (25, 30), (0.85, 1.15), (0.9, 1.1), (0.9, 1.1))