{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Fitting the data\nYou need to find a parameterized curve that can match all of these predictions. The simplest thing we can do is linear regression, and it seems to work fairly well with the known error of 70ml. We want to use a neural network to predict the slope and intercept for each FVC curve to do so we need to find the slope and intercept that best fits the data.\n\nTo deal with all of the issues associated with this fitting we can try to use bayesian programming as outlined in [this](https://docs.pymc.io/notebooks/GLM-robust.html) tutorial and as performed below. This could potentially allow us to find the bayesian estimate for confidence etc using a neural network,"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pydicom\nimport os\nimport random\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nfrom PIL import Image\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow as tf\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import pymc3 as pm\nimport theano","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def seed_everything(seed=2020):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    tf.random.set_seed(seed)\n    \nseed_everything(42)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT = \"../input/osic-pulmonary-fibrosis-progression\"\nBATCH_SIZE=128","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = pd.read_csv(f\"{ROOT}/train.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_ids = tr['Patient'].unique()\nfor pid in patient_ids[:20]:\n    mx = tr['Patient']==pid\n    plt.plot(tr['Weeks'].loc[mx],tr['FVC'].loc[mx])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_obs = []\nfor pid in patient_ids[:20]:\n    mx = tr['Patient']==pid\n    plt.plot(tr['Weeks'].loc[mx],tr['Percent'].loc[mx])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Fit the data and sample curves\nEach sample is scaled by a different amount in the below fit and so the slopes and intercepts are not comparable, they need to be transformed to the full range."},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_ids = tr['Patient'].unique()\nnsamples = 100\nnames = []; samples = np.zeros((2*len(patient_ids),nsamples))\nfor i in range(len( patient_ids )):\n# for i in range(3):\n    pid = patient_ids[i]\n    mx = tr['Patient']==pid\n    x_sample = tr['Weeks'].loc[mx]\n    y_sample = tr['FVC'].loc[mx]\n\n    xarr = np.array(x_sample)\n    xd = (xarr-min(xarr))/(max(xarr)-min(xarr))\n\n    data = dict(x=xd, y=np.array(y_sample))\n\n    with pm.Model() as model:\n        family = pm.glm.families.StudentT()\n        pm.glm.GLM.from_formula('y ~ x', data, family=family)\n        trace = pm.sample(nsamples, cores=2)\n    \n    plt.figure(figsize=(7, 5))\n    adjust = ((max(xarr)-min(xarr))+min(xarr))\n    mslp = trace['x'][-nsamples:]\n    aslp = mslp/adjust\n    mint = trace['Intercept'][-nsamples:]\n\n    plt.plot(x_sample, y_sample, 'x', label='data')\n    rng = np.linspace(min(xarr)-5,max(xarr)+5,11)\n    plt.plot(np.array([rng]*100).T,(np.outer(aslp[:100],rng)+mint[:100,np.newaxis]).T ,'k',alpha=0.05)\n    \n    samples[2*i:2*i+2] = np.vstack((aslp,mint))\n    names+=[pid]\n    \n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# import matplotlib as mpl\n# plt.hist2d(samples[0],samples[1])\n# plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.save('samples',samples)\nnp.save('names',np.array(names))","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}