Julia code in an example of XRjulia

Julia code in an example of XRjulia.ipynb

{
  "cells": [
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "See https://insightr.wordpress.com/2017/10/18/writing-julia-functions-in-r-with-examples/"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# = Define a function in Julia for a linear regression = #\n#regjl <- juliaEval(\"\nfunction reg(x,y)\n    n=size(x,1)\n    xreg=hcat(ones(size(x)[1],1),x)\n    k=size(xreg,2)\n    p1=((xreg'xreg)^(-1))\n    b=p1*xreg'y\n    r=y-xreg*b\n    sig=(r'r)/(n-k)\n    vmat=sig[1]*p1\n    sigb=sqrt.(diag(vmat))\n    t=b./sigb\n \n    return (b,t)\nend\n#\")\n \n# = Tell R regjl is a function = #\n#regjl_function=JuliaFunction(regjl)",
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 1,
          "data": {
            "text/plain": "reg (generic function with 1 method)"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# = Generate Data = #\n# set.seed(1)\n# x=matrix(rnorm(9000),300,30)\n# y=rnorm(300)\nsrand(1)\nx = randn(300,30)\ny = randn(300)\n\n# = Run Julia Regression = #\n#test=regjl_function(x,y)\n# = Convert Julia object to R = #\n#JLreg=juliaGet(test)[[1]]\n@time reg(x,y) \n\n# = Run regression using lm and compare = #\n#Rreg=coef(lm(y~x))\n#cbind(Rreg,JLreg)",
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": "  1.003914 seconds (571.21 k allocations: 28.702 MiB, 5.61% gc time)\n",
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "execution_count": 2,
          "data": {
            "text/plain": "([-0.0773863, -0.0313114, -0.0519369, 0.0753998, 0.0739596, 0.0263702, -0.000919084, 0.0231089, -0.0463216, -0.00319491  …  -0.0365581, 0.0870522, 0.034786, 0.0422321, 0.0232975, 0.0840548, -0.00566577, -0.041691, 0.0609428, 0.0289909], [-1.33399, -0.540728, -0.892494, 1.28071, 1.30626, 0.474433, -0.0163147, 0.413931, -0.857235, -0.0593057  …  -0.62097, 1.48257, 0.603138, 0.728385, 0.418197, 1.55388, -0.0989099, -0.751832, 0.998715, 0.473561])"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# = Load Combinatorics package in Julia. You need to install it. = #\n# = In julia, type: Pkg.add(\"Combinatorics\")\n#comb=juliaEval(\"using Combinatorics\")\nusing Combinatorics\n \n# = Evaluate the CSR function = #\n#csrjl <- juliaEval(\"\nfunction csr(x,y,k,K,fixed_controls)\n    fixed_controls=floor(Int,fixed_controls)\n    if fixed_controls!=0\n        w=x[:,fixed_controls]\n        nonw=setdiff(1:size(x,2),fixed_controls)\n    end\n    if fixed_controls[1]==0\n        nonw=1:size(x,2)\n    end\n\n    save_stat=zeros(size(x,2),2)\n    save_stat[:,2]=1:size(x,2)\n    for i in nonw\n        if fixed_controls[1]==0\n            (betas,tstat)=reg(x[:,i],y)\n            save_stat[i,1]=abs(tstat[2])\n        end\n        if fixed_controls!=0\n            (betas,tstat)=reg(hcat(w,x[:,i]),y)\n            save_stat[i,1]=abs(tstat[end])\n        end\n    end\n\n    t_ord = sortperm(save_stat[:,1],rev=true)\n    save_stat=save_stat[t_ord,:]\n\n    selected=save_stat[1:K,2]\n\n    #comb=collect(Base.combinations(selected,k))\n    comb=collect(combinations(selected,k))\n    m=size(comb)[1]\n    final_coef=zeros(m,size(x,2))\n    final_const=zeros(m)\n\n    if fixed_controls!=0\n        for i in 1:m\n            xr=hcat(x[:,fixed_controls],x[:,floor(Int,comb[i])])\n            (model,tstat)=reg(xr,y)\n            final_const[i]=model[1]\n            model1=model[2:end]\n            if length(fixed_controls)>1\n                final_coef[i,fixed_controls]=model1[1:(size(fixed_controls)[1])]\n                model2=model1[size(fixed_controls)[1]+1:end]\n            else\n                final_coef[i,fixed_controls]=model1[1]\n                model2=model1[2:end]\n            end\n\n            final_coef[i,floor(Int,comb[i])]=model2\n        end\n    elseif fixed_controls[1]==0\n        for i in 1:m\n            # (model,tstat)=reg(x[:,floor(Int,comb[i])],y)\n            (model,tstat)=reg(x[:,floor.(Int,comb[i])],y)\n            final_const[i]=model[1]\n            # final_coef[i,floor(Int,comb[i])]=model[2:end]\n            final_coef[i,floor.(Int,comb[i])]=model[2:end]\n        end\n    end\n    aux=hcat(final_const,final_coef)\n    result=[mean(aux[:,i]) for i in 1:size(aux)[2]]\n\nend\n#\")\n# = Define it as a function to R = #\n#csrjl_function=JuliaFunction(csrjl)",
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 3,
          "data": {
            "text/plain": "csr (generic function with 1 method)"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "@time csr(x,y,4,20,0);",
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": "  0.968043 seconds (668.74 k allocations: 194.576 MiB, 3.55% gc time)\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "ENV[\"LINES\"] = 100\n@time csr_Julia = csr(x,y,4,20,0)",
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "text": "  0.156353 seconds (184.67 k allocations: 169.663 MiB, 15.04% gc time)\n",
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "execution_count": 5,
          "data": {
            "text/plain": "31-element Array{Float64,1}:\n -0.0678484 \n -0.00691596\n -0.0095239 \n  0.0151644 \n  0.0131675 \n  0.0       \n  0.00348819\n  0.00683737\n -0.0112112 \n  0.0       \n  0.0       \n  0.00927163\n  0.0099952 \n -0.0119132 \n  0.0128721 \n  0.00751387\n  0.0       \n  0.0       \n  0.0       \n  0.00840607\n  0.0       \n -0.00746758\n  0.016471  \n  0.00923173\n  0.00923969\n  0.0       \n  0.019383  \n  0.0       \n -0.00863687\n  0.0142173 \n  0.0       "
          },
          "metadata": {}
        }
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "julia-0.6",
      "display_name": "Julia 0.6.0",
      "language": "julia"
    },
    "toc": {
      "threshold": 4,
      "number_sections": true,
      "toc_cell": false,
      "toc_window_display": false,
      "toc_section_display": "block",
      "sideBar": true,
      "navigate_menu": true,
      "moveMenuLeft": true,
      "widenNotebook": false,
      "colors": {
        "hover_highlight": "#DAA520",
        "selected_highlight": "#FFD700",
        "running_highlight": "#FF0000",
        "wrapper_background": "#FFFFFF",
        "sidebar_border": "#EEEEEE",
        "navigate_text": "#333333",
        "navigate_num": "#000000"
      },
      "nav_menu": {
        "width": "252px",
        "height": "12px"
      }
    },
    "language_info": {
      "file_extension": ".jl",
      "name": "julia",
      "mimetype": "application/julia",
      "version": "0.6.0"
    },
    "gist": {
      "id": "",
      "data": {
        "description": "Julia code in an example of XRjulia",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}