{"cells":[{"cell_type":"markdown","metadata":{"id":"_2gK8MYFDNHq"},"source":["Developed by Megan Renz, Rebeckah Fussell, and Natasha Holmes for Cornell Physics Labs."]},{"cell_type":"markdown","metadata":{"id":"Iv0-megcDNHs"},"source":["### Fitting"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"KuFrjSYbDNHs"},"outputs":[],"source":["import sys\n","\n","%matplotlib inline\n","#%matplotlib notebook\n","import matplotlib.pyplot as plt\n","import numpy as np\n","import csv\n","from scipy.optimize import minimize, rosen, rosen_der, curve_fit\n","from matplotlib.widgets import TextBox\n","import pandas as pd\n","\n","# def chiSquared2(args,x, y, dy, f ):\n","#     '''Function Chi-Squared.\n","#     x, y and dy are numpy arrays, referring to x, y and the uncertainty in y respectively.\n","#     f is the function we are fitting.\n","#     args are the arguments of the function we have fit.\n","#     '''\n","#     return 1/(len(x)*np.sum((f(x, args)-y)**2/dy**2))\n","\n","def chiSquared(x, y, dy, f, args):\n","    '''Function Chi-Squared.\n","    x, y and dy are numpy arrays, referring to x, y and the uncertainty in y respectively.\n","    f is the function we are fitting.\n","    args are the arguments of the function we have fit.\n","    '''\n","    return 1/(len(x))*np.sum((f(x, args)-y)**2/dy**2)\n","\n","def linear(x,args):\n","    '''\n","    A special case of Poly.\n","    '''\n","    return x*args[0]+args[1]\n","\n","def linearFit(x, m,b):\n","    return m*x+b\n","\n","def poly(x, args):\n","    '''\n","    returns the value of the polynomial sum (x**i*args[i])\n","    '''\n","    total=x**0*args[0]\n","    for i in range(1,len(args)):\n","        total+=x**i*args[i]\n","    return total\n","\n","def polyFit(x, A, B, C):\n","  return A*(x**2)+B*x+C\n","\n","def autoFit(x=[], y=[], dy=[], title=\"Use title= in your call\", xaxis=\"Use xaxis= in your call\", yaxis=\"Use yaxis= in your call\"):\n","    x0=[0,0]\n","    x=np.array(x)\n","    y=np.array(y)\n","    dy=np.array(dy)\n","\n","    args,cov=curve_fit(linearFit, x, y, x0, dy)\n","\n","    m=args[0]\n","    b=args[1]\n","    dm = np.sqrt(cov[0][0])\n","    db = np.sqrt(cov[1][1])\n","\n","    fig, ax=plt.subplots(1,2, figsize=(12,6), dpi=75)\n","    line,=ax[0].plot(x,linear(x, [m,b]))\n","    data=ax[0].errorbar(x,y, dy, fmt='.k', capsize=5, markersize=10)\n","    ax[0].set_title(title, fontsize=20)\n","    ax[0].set_xlabel(xaxis, fontsize=15)\n","    ax[0].set_ylabel(yaxis, fontsize=15)\n","    ax[1].set_title(\"Residuals\",fontsize=20)\n","    ax[1].set_xlabel(xaxis,fontsize=15)\n","    ax[1].set_ylabel(\"point -- line\",fontsize=15)\n","    residuals=y-linear(x, [m,b])\n","    res=ax[1].errorbar(x,residuals, dy, fmt='.k', capsize=5, markersize=10)\n","    ax[1].grid(True, which='both')\n","    plt.show()\n","\n","    print(\"chi-squared value =  {}\".format(chiSquared(x,y, dy, linear, [m, b]).round(3)))\n","    print(\"Best fit line :  \")\n","    print(\"y=\", round(args[0],4),\"x\", \"+\", round(args[1],4))\n","    print(\"m={}+\\\\-{}\".format(m.round(3),dm.round(3),))\n","    print(\"b={}+\\\\-{}\".format(b.round(3),db.round(3),))\n","\n","def manualFit(x=[], y=[], dy=[], m=[], b=[], title=\"Use title= in your call\", xaxis=\"Use xaxis= in your call\", yaxis=\"Use yaxis= in your call\"):\n","    x0=[0,0]\n","    x=np.array(x)\n","    y=np.array(y)\n","    dy=np.array(dy)\n","\n","    fig, ax=plt.subplots(1,2, figsize=(12,6), dpi=75)\n","    line,=ax[0].plot(x,linear(x, [m,b]))\n","    data=ax[0].errorbar(x,y, dy, fmt='.k')\n","    ax[0].set_title(title, fontsize=20)\n","    ax[0].set_xlabel(xaxis, fontsize=15)\n","    ax[0].set_ylabel(yaxis, fontsize=15)\n","    ax[1].set_title(\"Residuals\",fontsize=20)\n","    ax[1].set_xlabel(xaxis,fontsize=15)\n","    ax[1].set_ylabel(\"point -- line\",fontsize=15)\n","    residuals=y-linear(x, [m,b])\n","    res=ax[1].errorbar(x,residuals, dy, fmt='.k')\n","    ax[1].grid(True, which='both')\n","    plt.show()\n","\n","    print(\"chi-squared value =  {}\".format(chiSquared(x,y, dy, linear, [m, b]).round(3)))\n","    print(\"Fit line :  \")\n","    print(\"y=\", m,\"x\", \"+\", b)\n","\n","\n","def autoFit2(x=[], y=[], dy=[], title=\"Use title= in your call\", xaxis=\"Use xaxis= in your call\", yaxis=\"Use yaxis= in your call\"):\n","    x0=[0,0]\n","    x=np.array(x)\n","    y=np.array(y)\n","    dy=np.array(dy)\n","\n","    args,cov=curve_fit(linearFit, x**2, y, x0, dy)\n","\n","    m=args[0]\n","    b=args[1]\n","    dm = np.sqrt(cov[0][0])\n","    db = np.sqrt(cov[1][1])\n","\n","    fig, ax=plt.subplots(1,2, figsize=(12,6), dpi=75)\n","    line,=ax[0].plot(x,linear(x**2, [m,b]))\n","    data=ax[0].errorbar(x,y, dy, fmt='.k')\n","    ax[0].set_title(title, fontsize=20)\n","    ax[0].set_xlabel(xaxis, fontsize=15)\n","    ax[0].set_ylabel(yaxis, fontsize=15)\n","    ax[1].set_title(\"Residuals\",fontsize=20)\n","    ax[1].set_xlabel(xaxis,fontsize=15)\n","    ax[1].set_ylabel(\"point -- line\",fontsize=15)\n","    residuals=y-linear(x**2, [m,b])\n","    res=ax[1].errorbar(x,residuals, dy, fmt='.k')\n","    ax[1].grid(True, which='both')\n","    plt.show()\n","\n","    print(\"chi-squared value =  {}\".format(chiSquared(x**2,y, dy, linear, [m, b]).round(3)))\n","    print(\"Best fit line :  \")\n","    print(\"y=\", round(args[0],4),\"x^2\", \"+\", round(args[1],4))\n","    print(\"m={}+\\\\-{}\".format(m.round(3),dm.round(3),))\n","    print(\"b={}+\\\\-{}\".format(b.round(3),db.round(3),))\n"]},{"cell_type":"markdown","metadata":{"id":"kxQYMbaqDNHt"},"source":["## Uncertainty and $t^\\prime$"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"OTHBWlSTDNHu"},"outputs":[],"source":["def standard_deviation(data):\n","    '''\n","    Takes in a one-dimensional numpy array argument, and\n","    returns the standard deviation of the dataset.\n","    '''\n","    N = len(data)\n","    return np.sqrt(np.sum((data - np.mean(data))**2)/(N-1))\n","\n","def standard_unc_of_mean(data):\n","    '''\n","    Takes in a one-dimensional numpy array argument, and\n","    returns the standard uncertainty of the mean of the dataset.\n","    '''\n","    N = len(data)\n","    return standard_deviation(data) / np.sqrt(N)\n","def t_prime(A, dA, B, dB):\n","    '''\n","    Returns the value of t prime. Requires 4 arguments in this order:\n","    measurement A, uncertainty of measurement A, measurement B, uncertainty of measurement B.\n","    '''\n","    return abs(A - B) / np.sqrt(dA**2 + dB**2)\n","\n","\n","\n","\n"]},{"cell_type":"markdown","metadata":{"id":"6R-8HmsODNHu"},"source":["### Uncertainty propagation"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"gJR-BpG2DNHu"},"outputs":[],"source":["def dsum(x, dx, y, dy):\n","    return np.sqrt(dx**2+dy**2)\n","\n","def dsub(x, dx, y, dy):\n","    return np.sqrt(dx**2+dy**2)\n","\n","def dmult(x, dx, y, dy):\n","    return x*y*np.sqrt((dx/x)**2+(dy/y)**2)\n","\n","def ddiv(x, dx, y, dy):\n","    return x/y*np.sqrt((dx/x)**2+(dy/y)**2)\n","\n","def dpwr(x, dx, n):\n","    return n*dx*(x**n)/x\n","\n","def dln(x, dx):\n","    return dx/x"]}],"metadata":{"kernelspec":{"display_name":"Python 3 (ipykernel)","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.9.13"},"colab":{"provenance":[]}},"nbformat":4,"nbformat_minor":0}