Error calculation with linear regression

You will want to first convert your current uncertainties $\delta i_c$ into uncertainties in $\log{i_c}$, $\delta(\log{i_c})$. You will need the uncertainties in the $x$ and $y$ values that you use in the fit in order to get the uncertainties in the fit parameters. Note that the uncertainty of the $\log$is not simply the $\log$ of the uncertainties.

Taylor derives the formulas for the uncertainties in the fit parameters in the case that all of the underlying measurements have the same uncertainties. He leaves the case where the uncertainties of the underlying measurements are not the same to the problems for the reader to derive. The derivation is not difficult for either case, but it is tedious for both of them. I don't think it's problematic to simply quote the results for you; they are well-known, and I'm assuming that you're just supposed to look them up. Frankly, I think they should be given to you by your fitting software, but everything I've used makes it difficult to find them.

Using your notation for the slope $\alpha$ and the y-intercept $b$, the uncertainties are: $$ \sigma_\alpha = \sqrt{\frac{\sum{wx^2}}{\Delta}}, \sigma_b = \sqrt{\frac{\sum{w}}{\Delta}} $$ where $x$ is the set of your $V_{BE}$'s, $w = 1/(\sigma_y)^2$, $\sigma_y = \delta(\log{i_c})$, $\Delta = N\sum{x^2} - (\sum{x})^2$, $N$ is the number of points you are using for the fit, and the summations are over all the points you are using for the fit.

When you have the uncertainty for each of your four values of $k_b$, you can get the uncertainty for their mean by adding the four uncertainties in quadrature, like you would for the uncertainty of a sum of measurements: $\sigma_\mathrm{total} = \frac{1}{4}\sqrt{\sigma_1^2 + \sigma_2^2 + \sigma_3^2 + \sigma_4^2}$.

For a linear regression fit of the form: $y=ax+b$, with $S_{yi}=S_i\neq const.$, then: $$ a=\frac{\sum{\frac{1}{S_i^2}} \sum{\frac{x_iy_i}{S_i^2}} - \sum{\frac{x_i}{S_i^2}} \sum{\frac{y_i}{S_i^2}}} {\sum{\frac{1}{S_i^2}} \sum{\frac{x_i^2}{S_i^2}} - \left(\sum{\frac{x_i}{S_i^2}}\right)^2} $$ with: $$ S_a=\sqrt{\frac{\sum{\frac{1}{S_i^2}}} {\sum{\frac{1}{S_i^2}}\sum{\frac{x_i^2}{S_i^2}}-\left(\sum{\frac{x_i}{S_i^2}}\right)^2}} $$ And from your expression for $\alpha$ you can determine the $S_\alpha$.

If your data has error bars for the $y$-axis, you could use a least-squares/$\chi^2$ technique to fit your model:

You have $\{x_i,y_i\pm\epsilon_i\}$ for $i=1,\ldots,N$ data. You want to find $K_B$ and $I_c$ such that $y=f(k_B,I_C,x)$. Construct the function, $$ \chi^2 = \sum_{i=1}^N\frac{(f(k_B,I_c,x_i)-y_i)^2}{\epsilon_i^2}, $$ and minimize this function wrt to $k_B$ and $I_s$ to find your best estimates for these parameters, say $k_B^*$ and $I_s^*$. To find the confidence intervals for your estimates (at 68%), fix $I_s=I_s^*$, and vary $k_B$ up and down until $\delta \chi^2 = \chi^2(k_B,I_c^*) - \chi^2(k_B^*,I_c^*) = 1$. The interval for which $\delta \chi^2 \le 1$ is the 68% confidence interval for $k_B$.

This works because Wilkes' theorem shows that $\delta \chi^2$ is asymptotically $\chi^2$-distributed.